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ABSTRACT 

In recent years, there has been a proliferation of wide-field sky surveys to search for a variety of 
transient objects. Using relatively short focal lengths, the optics of these systems produce undersam- 
pled stellar images often marred by a variety of aberrations. As participants in such activities, we 
have developed a new algorithm for image subtraction that no longer requires high quality reference 
images for comparison. The computational efficiency is comparable with similar procedures currently 
in use. The general technique is cross-convolution: two convolution kernels are generated to make a 
test image and a reference image separately transform to match as closely as possible. In analogy to 
the optimization technique for generating smoothing splines, the inclusion of an RMS width penalty 
term constrains the diffusion of stellar images. In addition, by evaluating the convolution kernels 
on uniformly spaced subimages across the total area, these routines can accomodate point spread 
functions that vary considerably across the focal plane. 

Subject headings: techniques: photometric, methods: statistical 
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1. INTRODUCTION 

The advent of low-noise megapixel electronic image 
sensors, cheap fast computers, and terabyte data storage 
systems has enabled searches for rare astrophysical phe- 
nomena that would otherwise be effectively undetectable. 
Examples include the discoveries of MACHOs, transits of 
extra-solar planets, and vastly greater numbers of super- 
novae. Common to all of these efforts is a need to repeat- 
edly image large portions of the sky to uncover rare and 
subtle changes of brightness. This forces the observer to 
contend with images taken under less than ideal condi- 
tions such as poor weather and crowded star fields. To 
solve the problem of comparing images taken with differ- 
ent seeing conditions, a number of research groups have 
developed image subtraction algorithms which compen- 
sate for image blurring effects prior to image differencing. 
In this paper, we describe a variation of this technique 
which treats pairs of images in a symmetric fashion, re- 
ducing the requirements of first obtaining an ideal refer- 
ence image. The software code is relatively simple and 
has been made freely available. 

Following our initial discovery of prompt optical radia- 
tion from a gamma-ray burst in 1999, our ROTSE collab- 
oration set out to construct a set o f four identical wide- 
field telescopes (jAkerlof et al.lf2003h to explore these phe- 
nomena more deeply. The resulting instruments, called 
ROTSE-III, were installed in Australia, Texas, Namibia 
and Turkey in 2003 and 2004. By explicitly choosing a 
fast (//1.9) optical system and short focal length (850 
mm), we provided the option to search for possible or- 
phan GRB afterglows without unduly compromising the 
potential for mapping GRB optical light curves at early 
times. Although we always intended to use these instru- 
ments for generic astrophysical optical transient searches, 
it was recognized that the plate scale (3.3" per pixel) was 
too coarse for easy identification of a supernova embed- 
ded in a normal galaxy. 
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This challenge was addressed by Robert Quimby when 
he became a graduate student at the University of Texas. 
As an undergrad at the University of California at Berke- 
ley, he had worked extensively with the Supernova Cos- 
mology Project (SCP) and was quite familiar with the SN 
discovery pro cess. Quimby adapted t he SCP image sub- 
traction code (jPerlmutter et all l999) for use as the basic 
tool for finding 30 SNe over a period of two years from 
observations with the University of Texas 30% alloca- 
tion of ROTSE-IIIb time at the McDonal d Observatory 
Amo ng those discoveri es are SN 2005apjQuimbv et al.1 
120071 ) and SN 2006gy ([Smith et alJl2007h which appear 
to be the intrinsically brightest SNe ever identified. 

In view of the e vident succes s of the Texas Super- 
nova Search (TSS) (|Quimbvll2006l ). our group at the Uni- 
versity of Michigan probed the image subtraction prob- 
lem with the goal of applying this to the considerably 
more extensive image data available to the entire suite 
of ROTSE-III telescopes. The original hope of using the 
SCP code was abandoned following the realization that 
the program would not be made freely available. We at- 
tempted to adopt the ISIS image subtraction package 2 
but were discouraged by the initial results. The signifi- 
cant undersampling of ROTSE-III stellar images coupled 
with asymmetric point spread functions across the image 
plane created a severe challenge for making clean sub- 
tractions. These issues are not satisfactorily ad d ressed 
by t he alg o rithm s described by lAlard fc Luptonl (|1998f ) 
andlAJardJ (|2000l) for two reasons: (1) we do not always 
have the luxury of a substantially higher quality refer- 
ence image and (2), the point spread functions (PSFs) 
are often approximately elliptical with the axes oriented 
at any angle in the image plane. For a variety of rea- 
sons, the ROTSE-III PSFs can vary with temperature 
and telescope orientation. Thus, the possibility that one 
can simply convolve a new image to an ideal reference im- 
age is not always viable. With this in mind, we sought 
to develop a more symmetric algorithm that would be 
robust enough to handle less pristine observations. We 

2 http:/ /www2. iap.fr/users/alard/package. html 
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should emphasize that the aim of this project is primar- 
ily for the reliable identification of transients in a very 
large database, not precision photometry. 

2. MATHEMATICAL METHOD 

The basic technique for image subtraction presented 
by Alard and Lupton depends on finding a suitable PSF 
smearing kernel, K{u, v) that when convoluted with the 
reference image, R(x,y), generates a transformed image, 
R*(x, y), that can be compared on a pixel- by-pixel basis 
with a new test image of lesser quality, T(x, y): 

R(x,y)®K(u,v) = R*(x,y) (1) 

The kernel, K(u,v), is constructed by a linear superpo- 
sition of basis functions of the form: 

fn,p, q ^^) = ^ q e~ {u2+v2)/2al (2) 

Alard and Lupton recommend a three-fold ensemble of 
terms with a n values spanning a 9-fold range. These 
functions are poor choices to synthesize an elliptical PSF 
at an arbitrary angle with respect to the imager sensor 
axes although they will be satisfactory for PSFs with 
close to azimuthal symmetry. The specific values for o~ n , 
p and q must be determined by ad hoc comparisons with 
the characteristic PSFs associated with the particular 
instrument in use. The set of amplitudes for the basis 
functions is computed by the least squares technique to 
minimize the pixel-by-pixel differences between R* and 
T. 

From this starting point, we decided to symmetrize 
the Alard-Lupton procedure by creating two convolution 
operators so that: 

R(x, y) <g> K R (u, v) w T(x, y) ® K T {u, v) (3) 

In the limit that the reference image is substantially bet- 
ter than the test image, the Kr operator smears R with 
the point spread function characteristic of T while K~t 
will be essentially equal to the identity operator so that 
its effect on T will be negligible. Under these conditions, 
the computation becomes functionally equivalent to the 
procedure adopted by Alard and Lupton. However, in 
general, the addition of a second convolution operator 
injects new mathematical degrees of freedom that must 
be constrained. The most obvious is that Kr and Kt 
can be multiplied by an arbitrary constant without vi- 
olating image convolution equality. This can be conve- 
niently resolved by demanding that one or both kernels 
be flux-conserving, ie: 

J2K(u,v) = l (4) 

The second, more complex, problem arises from the dif- 
fusion of a stellar image if Kr and Kt approximate broad 
Gaussian distributions. The convoluted image equality 
will be maintained but the signal-to-noise ratio of the 
subtracted image will drastically diminish. The solution 
to this was found by analogy to a similar problem in 
the application of smoothing splines to drawing curves 
through data with errors. In the latter case, one could 
trivially create a spline curve that ran exactly through 
each and every data point (as long as the abscissas are 
distinct). Such a curve would appear very wiggily and 



poorly represent the trend of the data. The solution to 
this problem is to add a curvature penalty term to the 
least squares residuals so that a trade-off is reached be- 
tween adequately fitting the data and inserting unneces- 
sarily complex behavior into the smoothed interpolation. 
The coefficient that scales the curvature term is a mea- 
sure of the stiffness of the spline. There exists an elegant 
method called cross-validation that determines the stiff- 
ness parameter from the standard deviation errors for 
each data point. 

For image subtraction, the degradation of the signal- 
to-noise ratio is proportional to the effective number of 
pixels that are summed by the convolution kernel. Since 
each pixel has an associated variance, a pix , the variance 
for a signal diffused over N pix pixels will be N pix a pix . 
Assuming Gaussian distributions, we can estimate N p i X 
from the width of the effective stellar point spread func- 
tion, N p i X ~ 4:Tr(ap SF + o- 2 K ), where apsF is the ba- 
sic stellar PSF width (in pixels) and ok is the diffusive 
width of the convolution, K. If <j 2 K is evaluated by the 
normal formula, o 2 K — K^r 2 where Ki are the kernel 
element amplitudes, the sum can be deceptively small 
when the values for Ki alternate in sign. Noting that 
(Ki) cx Jj-, the value for a 2 K can be better estimated 

by 51F Kfrf which equally penalizes both positive and 
negative contributions to the kernel elements. Although 
the scaling behavior of the penalty coefficient is under- 
stood in terms of the image size, o~ pix , and crp SF , we 
have not investigated whether there is an elegant way to 
evaluate this quantity analogously to the cross-validation 
technique for splines. 

3. COMPUTATIONAL METHODS 

The image preprocessing that we require is similar to 
that described by Alard and Lupton. Flat-fielded ima ges 
are processed by SExtractor (jBertin &: Arnoutil 19961 ) to 
create object lists with precise stellar coordinates. The 
IDL routine 3 , POLYWARP, is used to warp the new test 
image to overlay stellar objects in the reference image 
as closely as possible. A valid pixel map is generated 
to avoid pixels close to the image perimeter and screen 
against saturated values, etc. At this point, the fun- 
damental image subtraction code is invoked as an IDL 
routine which first performs image flux normalization to 
equalize the mean values of the two images under com- 
parison. 

The most basic choice that the user must make is the 
representation of the convolution kernels. We have re- 
stricted them to n x n arrays with n odd. This permits a 
simple representation for the convolution identity opera- 
tor: -ftT[„/2].[n/2] — 1 while all other elements of Kij are 
zero. For the ROTSE-III images, n = 9 appears to pro- 
vide more than adequate coverage of stellar point spread 
functions under the worst conditions. 

The values for the convolution kernel elements are de- 
rived from the difference image: 

D(x, y) = {R(x, y) <g> K R (u, v)} - {T(x, y) <gs K T (u, »)} 

.(5) 

Invoking the criterion that D{x, y) should be a mini- 
mum subject to the requirements for Kr and Kt gener- 

3 ITT Visual Information Solutions, ITT Industries, Inc. 
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ates 2(n 2 — 1) linear equations via the usual least squares 
procedure to solve for the independent coefficients for 
Kr(u, v) and Kt(u, v) after imposing the kernel unitar- 
ity constraints. As described earlier, these equations will 
not provide unique solutions for Kr and Kt because the 
effective width of the two convolution kernels can still 
be radially scaled without substantially affecting the dif- 
ference image, D. Thus, the quantity to be minimized 
must include a penalty term for radially diffusing the 
convoluted images any further than necessary. Following 
from earlier remarks, this figure-of-merit function can be 
represented as: 

Q = Y1 D ( x > y) 2 +xY,(u 2 +v 2 ) 2 {K R (u, vf+K T {u, v) 2 } 

(6) 

where A is a constant selected to balance the contribu- 
tions of the two competing error terms. From the discus- 
sion given above, the value for A should scale as: 

A = 2irN lmage - R 2 — A' (7) 

a PSF 

where Ni mage is the total number of pixels in the image, 
a 2 R and a 2 -, are the pixel amplitude variances, opsf is 
the characteristic stellar PSF width and A' is a constant 
of order unity. 

With two 9x9 convolution kernels, the number of free 
parameters is 160 and the size of the regression matrix 
becomes problematic. The main concern is that if the 
images are essentially featureless (ie. no stars), the ma- 
trix elements become indistinguishable and the inverse 
matrix will be singular. To avoid these effects as well as 
various other computational issues, a binary valued mask 
array is created to eliminate sampling around the image 
perimeter, saturated pixels and all featureless areas not 
associated with stellar objects as determined by SExtrac- 
tor. This approach was quite successful: the degree of 
singularity of the regression matrix was determined dur- 
ing the inversion process using the IDL singular value 
decomposition routines SVDC and SVSOL code s derived 
from Numerical Recipes in C ()Press et al.lll992l ). 

For our ROTSE project, computational efficiency is 
critical because we typically acquire 400 images per night 
with each telescope and these must be reduced in situ 
comfortably within 24 hours. It was easily verified that 
most of the image subtraction calculations were devoted 
to computing the convolution kernels regression matrix 
described above. Examination of the two-dimensional 
structure of the kernels showed that the amplitudes near 
the edges of the 9x9 arrays were always small and sug- 
gested that the representation could be significantly re- 
duced from 81 values to 25 by assuming a mapping from 
a reduced number of wavelet functions. Thus, each con- 
volution kernel was represented by a linear superposition: 

K(u,v) =^2A. 1 B z (u,v) 

with the basis functions, Bi, chosen as discrete approx- 
imations to bicubic spline functions with characteristic 
widths of 1, | and 2 pixels, centered as shown in Fig- 
ure 1. Using this technique shrank the regression matrix 
from 160 x 160 to 48 x 48 with a consequent reduction 
in processing time of about an order of magnitude. This 
brought the computation throughput to values similar to 



what Robert Quimby had obtained using the SCP code 
as adapted for his Texas Supernova Search. 

Most of the image subtraction code was written in IDL 
with the exception of the evaluation of the regression ma- 
trix. Since this is the core of the computational burden 
and IDL is not particularly efficient in handling the nec- 
essary array indexing, this portion was coded in C and 
linked to the rest of the IDL programs using the IDL 
standard external calling interface. A crucial detail, par- 
ticularly important for the ROTSE-III telescopes, is the 
variation of the stellar PSF across the image plane. To 
accomodate this problem, each 2045 x 2049 image was 
subdivided into 36 sub-images of roughly equal size. The 
set of 50 kernel amplitude coefficients was calculated, one 
by one, for each of these sub-images. Cross-convolved 
images were obtained by bilinear interpolation for ev- 
ery pixel using the four nearest neighbor coefficient sets. 
The unitarity of the convolution kernels is guaranteed by 
the linearity of the interpolation method with respect to 
the gridded coefficient reference values. Although this 
sounds somewhat complicated, the calculation was ex- 
tremely efficient. 

Another significant concern is the estimation of the 
background sky intensity. Initially, we relied on SExtrac- 
tor to remove this background before subtraction. How- 
ever, when applying our code to images containing large 
galaxies such as M31 and M33, we realized that these 
backgrounds are poorly estimated around the cores of 
bright galaxies. The solution we adopted only removes 
the background difference between the two images in- 
stead of the individual background for each image sepa- 
rately. After the images are scaled so that stellar fluxes 
match, a sky difference map is generated by first perform- 
ing pixel by pixel subtraction. The low frequency spatial 
variation of this difference image is obtained by a process 
similar to the one used by SExtractor. The difference im- 
age is divided into 32 x 32 pixel subimages and median 
pixel values are recursively evaluated, subject to the con- 
straint that pixels with 3er excursions from the median 
are ignored. Saturated pixels are also excluded from the 
median computation. The resultant slowly varying back- 
ground is subtracted from one of the input images before 
invoking the cross-convolution algorithm. The remaining 
common non-zero background doesn't affect the estima- 
tion of the convolution kernels and the final subtracted 
image will be background-free. 

A comparison of the results of the cross-convolution 
and the single convolution algorithms is shown in Fig- 
ure 2. In the limiting case where the PSFs of the refer- 
ence image are azimuthally symmetric, the two meth- 
ods should produce rather similar results. However, 
when that condition is not satisfied, the cross-convolution 
method is more appropriate. 

For anyone wishing to employ the cross- 
convolution technique described in this paper, the 
source code can be dow nloaded from the URL, 
|http://hdl.handle.net/2027.42/57484[ which points to 
a set of files stored in Deep Blue, the University of 
Michigan institutional repository. Note: This Web 
site is currently inaccessible. Please contact the 
authors directly. 

4. OPERATIONAL EXPERIENCE 
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Subtraction of a 2045 x 2049 ROTSE image from a 
reference frame using the method described above takes 
approximately 4 minutes with a 2.0 GHz personal com- 
puter. If the same reference image is used multiple times, 
it needs to be convolved with the base kernels just once, 
saving computational time. Subtractions of three typical 
ROTSE images from the same reference frame takes ~ 
10 minutes on the same processor. It should be noted 
that the memory allocation for the process, mainly for 
storing the base kernel convoluted images, scales with the 
size of the image and number of kernel basis sets. For 
our choice of 25 kernel base vectors and a 2045 x 2049 
image size, ~ 1 GB memory is required. This is not a 
serious handicap for a modern desktop computer. 

Since August 2007, a supernova search pipeline us- 
ing this subtraction code has been running on images 
taken by the ROTSE-IIIb telescope. Selected fields with 
nearby rich clusters and a high density of known galax- 
ies are monitored on a daily basis, weather permitting, 
to a typical limiting magnitude of 18.5. For each field, 
two sets of four 60-sec exposures (20-sec for fields with 
bright target galaxies) are taken with a 30-minute ca- 
dence. Following the method developed by the Texas 
Supernova Search, the images for each 4-exposure epoch 
are co- added as are the total 8 images for the night. All 
three co-additions are subtracted by the same reference 
image. The difference images are processed through SEx- 
tractor to find residual objects. To reject false detections 
due to bad pixels, cosmic rays, asteroids and subtraction 
noise, further filtering is applied. The signal-to-noise ra- 
tio of a candidate has to be above 5 in the nightly 8-fold 
sum and 2.5 for the sum of single epoch. The positions of 
a candidate in 3 subtractions must match to within one 
pixel for detections above 15 SNR and 1.5 pixels for those 
with SNR below 15. The FWHM of the candidate has 
to lie within the range of one pixel and twice the median 
FWHM for stars in the convolved reference image. Fi- 
nally, minimum flux change cuts are applied with a lower 
threshold for detections embedded in known galaxies and 



higher for those corresponding to stellar objects. This 
later criterion is intended to suppress variable stars. 

In the 5-month period to date, the pipeline has iden- 
tified all 13 reported supernovae that lie within our 
searched fields. One of these initially escaped but was 
detected following modification of the mask size to pro- 
vide better performance during bad seeing. Also due 
to our early inexperience, two of these SNe in relatively 
bright galaxies were initially missed during hand scan- 
ning. In addition, the pipeline detected 7 novae in the 
fields of M31 and M33. Two novae rather close to the 
center of M31 were missed before the background eval- 
uation problem was addressed as described in section 3. 
In terms of transient recovery efficiency, both real- world 
and limited Monte Carlo comparisons show that our sub- 
traction code is comparable to the modified version of the 
Supernova Cosmology Project search code employed by 
the Texas Supernova Search. 

5. SUMMARY 

The algorithm described in this paper can be adapted 
for a wide variety of photometric searches for transient 
objects. Its performance appears to be at least as good 
as other codes presently in use. Since the method is de- 
signed to handle images with significantly varying qual- 
ity, it should remain effective when alternative programs 
may fail. 
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Fig. 1. — Diagram of the location of the 25 bicubic B-splines used to construct the convolution kernels. The 9 circles, 8 squares and 8 
hexagons mark the centers of the B-spline maxima with widths of 1, | and 2 pixels respectively. The dashed lines indicate the 9x9 grid 
of the underlying convolution kernels. 
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Fig. 2. — Comparison of image subtractions using the cross-convolution method described in this paper and the single convolution 
method described by Alard and Lupton and implemented in the ISIS code. The initial images were obtained by the ROTSE-IIIb telescope 
at McDonald Observatory. Shown here are 260 X 260 pixel subframes centered on a = 16 : 50 : 02.21, <5 = +23 : 46 : 32.88, covering a field 
of 0.235° X 0.235°. To demonstrate the results, three artificial "variable" stars were added to the test image (a) and the reference image (b) 
with PSFs appropriately matched to their respective fields. The locations are shown by black arrows. The subtracted image obtained by 
cross-convolution is depicted in (c) and the Alard-Lupton results arc shown in (d). The bright star near the lower right corner of the images 
has been replaced with a uniform gray level since neither subtraction technique can extract useful information from saturated pixels. 



